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Abstract. We consider the approximation of a class pseudodifferential operators by sequences of op- 
erators which can be expressed as compositions of differential operators and their inverses. We show that 
t e error m such approximations can be bounded in terms of the L x error in approximating a convolution 
kernel, md use this fact to develop convergence results. Our main result is a finite time convergence anal- 
ysis of the Engquist-Majda [7] Pade approximants to the square root of the d’Alembertian. We also show 
that no spatially local approximation to this operator can be convergent uniformly in time. We propose 
some temporally local but spatially nonlocal operators with better long time behavior. These are based on 
Itaguerre and exponential series. 
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1. Introduction. Pseudodifferential operators in time and/or space arise naturally in 
the formulation of many wave propagation problems. Prime examples include the restriction 
of unbounded domains via the introduction of artificial boundaries [7], the reduction of 
scattering problems to boundary integral equations [13], the factorization of wave operators 
into “one-way” equations [2], and the modeling of systems with memory [14]. Although there 
have been advances in recent years in the efficient application of certain integral operators, 
a typical approach to the numerical approximation of these problems is to replace the' 
pseudodifferential operator by a combination of differential operators and their inverses, 
which can then be applied using standard numerical techniques. This is the approach 
discussed in all the works cited above. 

In this note we examine the convergence of sequences of local approximations to a class 
of pseudodifferential operators. Our main example is the square root of the d’Alembertian, 
which is important in the study of acoustic, elastic and electromagnetic waves. Although 
many approximations to this operator have been proposed (e.g. [19]), error estimates of 
the type considered here do not seem to have been derived. (For estimates based on high- 
frequency asymptotics see [12].) Our estimates are important if the use of high order 
conditions is considered, as in [5]. We express the operator as the sum of a differential 
operator and a convolution operator with an L x kernel. Decomposing the local approxima- 
tion in the same way, we estimate the error in terms of the difference between the exact 
and approximate convolution kernels. We see that a sequence of local approximations is 
convergent if its sequence of approximate kernels converges (in L ± ) to the exact kernel. 

The main result we obtain in this framework is a finite time convergence theorem for 
the Pade approximants to the square root of the d’Alembertian which were proposed by 
Engquist and Majda [7]. We also show that this approximation (or any other which is 
local in space as well as time) is not convergent in our sense for infinite time. We propose 
some new approximations for infinite time which are spatially nonlocal. They are based on 
Laguerre and exponential series. Although we have not proven convergence for these, we 
do show that the error is reasonably small uniformly in time. The practical implementation 
and testing of long-time approximations will be carried out elsewhere. 

2. A Class of Pseudodifferential Operators and Local Approximations. We 
consider first an operator, A, acting on functions of t e [0, oo) which vanish along with their 
derivatives at t = 0. We assume: 

Assumption 1. 

. ^ d^u 

= €^([ 0 , 00 )). 


Making a Laplace transformation we then obtain: 

( 2 - 0A ) A(s) = a jS j + a_i(s), 

i— o 

where a_i (s) is the Laplace transform of the function (t). Typically, we know A directly 

To verify that Assumption 1 holds we must check that the inverse transform of a_ x is in 
Li. Necessary and sufficient conditions for this to hold are given in [20, Ch. 7]. 
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We give two examples, the first of which is relevant to problems involving the wave 
equation and the second of which is relevant to problems involving the diffusion equation. 
Example: Let 


(2.0.2) .4(3) — y / S 2 _J_ 1 ~ 5 y-- ■ 

5 -f \/ s 2 -f 1 

Then [15] 

(2.0.3) Au=~+S*u, S(t) = 

at t 

where Ji(t) is a first order Bessel function. It is easily verified that S 6 Li(f0 00 )) so 
Assumption 1 holds. 

Example: Let 

( 2 -0.4) A(s) = ^/7+T. 

Then A does not satisfy Assumption 1. 

We remark that the second example might be studied by expressing A as the compo- 
sition of a differential operator and a convolution. For example, 

( 2 - 0 - 5 ) Vs + 1 = A — ( s + !)> 

\/S+T ’ 

implies 

(!•...) *.(,.(*4.). TO.iJ. 

2.1. Local Approximations. We now consider local approximations B to operators 
A satisfying Assumption 1: 


(2.1.7) 


Bu ~ a 


i = 0 


dPu 

dP 


-f -S— x * u, 


where we assume that the differential operators in our standard decomposition of A and B 
agree. So that B will be local, we restrict B_x(t) to the class 71 defined by: 

Definition 1. A function f is in class 7Z if its Laplace transform, f(s), is a rational 
function. 


Alternatively, 71 may be described through the following elementary result. (See [15].) 

Lemma 2.1. The class 71 consists of functions which are products of polynomials, 
exponential functions and trigonometric functions. 

The error in approximating A by B is equivalent to the error in approximating convo- 
lution by .A_i by convolution by £J j . We therefore have; 

Theorem 1. Let u e L p ([0, T}), 0 < T < oo. Then 
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Proof: Clearly 

( 2 * 1 * 8 ) IK- 4 ~ b Ml p ([0,t]) = IK-4-1 ~B- 1) * «|U J> ([o,T])- 

The result then follows from the basic inequality ||/ * g\\ Li> < f|/||i 1 ||p||i, .o 

By this result we see that the problem of constructing accurate/convergent local ap- 
proximations is equivalent to approximating the L 1 function by elements of TZ. We 
therefore have: 

Corollary 1 . Let 0 < T < oo. For any e > 0 there exists a local operator B such 
that: 


II- 4 - #IU p ([0,T]) < € - 

(Here we use the operator norm.) 

Proof: We need only show that 7 1 is dense in L lm For finite T this follows from the density 
of the polynomials. For T = oo we can use the density of certain exponential families, which 
follows from logarithmic mappings of [0, oo) to finite intervals.o 

In subsequent sections we will study the error for specific sequences of approximations. 
It is nonetheless interesting to consider the possibility of optimal approximations where 
the degree of the transform of B - 1 is fixed. For finite intervals and the more restricted 
class of functions defined by products of polynomials and exponentials, the existence of 
best approximations in L\ has been established [4, Ch. 6]. Moreover, for monotone kernels 
interpolants can be computed and the degree of approximation estimated. In [10], boundary 
conditions based on time interpolation of were proposed and shown to have good 

Li approximation properties. 

2.2. Operators in Higher Dimension. In most applications the operators to be 
approximated act on functions defined on the product of a time interval with a spatial 
region, Cl. In the simplest case we take fi to be a hyperplane or torus and assume that A 
is a homogeneous function of s and |fc| = (^- Ar?) 1 / 2 , where the kj 9 s are Fourier variables. 
For example, if A = j s the square root of the d’Alembertian we have: 

(2.2.9) 6V2 = J s * + |jfe|2 = 5 + 1 *1 z = — . 

2 + \/z 2 + 1 |&| 

The [fc (-dependence of the operator complicates the approximation in a number of 
ways. First of a 11, an approximation is local in space only if the coefficients of the rational 
transform involve only even powers of \k\. In the example above, the rational function of z 
whose inverse transform approximates S must have numerator and denominator of definite 
and opposite parity in order to be spatially local. 

Second, since |fc| is not bounded, the approximations must be estimated uniformly in |jfe|. 
It is even possible that approximations can lead to ill-posed problems for some applications. 
In [18] the well-posedness of local approximations to the square root of the d’Alembertian is 
thoroughly studied. For general (spatially nonlocal) approximations such complete results 
are unknown, but various special cases have been considered [6, 11]. 
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3. The Engquist- Majda_Pa.de Approximants. The best known and most used 
rational approximations to y/s 2 + 1 are the Fade approximants suggested by Engquist and 
Majda [7J. In our language they may be written as (see [19]); 


(3.0.10) 


1 . U n „i(is) _ . sinnfl 

5 + y/l + s 2 U n (is) %in(n + l)0’ 


cos 6 = is , 


where U n denotes the Chebyshev polynomial of the second kind. Here we give a completely 
different derivation of this approximation. We begin with the formula [1, Ch. 9]: 


(3.0.11) S — — ~ — = — f y/l — w 2 cos wt dw. 

t vj - 1 

For fixed t we appr oximate the integral by the Gaussian quadrature rule appropriate for 
the weight v 1 - w 2 , based on second kind Chebyshev polynomials [1, Ch. 25]: 


(3.0.12) 


Mt) 


fcn(t) = sin 2 f -'j cos /(cos — — -)t^j 

n + 1 ti Vn + 1/ V n + l J ) 

t 2n cos£t, -1<£<1. 


(2n!)2 2n + 1 

Laplace transformation in time leads to the formula: 


(3.0.13) 

Remarkably we have: 

Theorem 2. Approximations (3.0.10) and (3.0.13) are the same. 

Proof: We first note that both expressions represent rational functions with denominators 
of exact degree n and the same poles, s~icosg i ,l=l... n . Moreover, the numerators 
are of degree n- 1. To prove their equality we will multiply each expression by (s - i cos ) 

and take the limit s^i cos ^ . For (3.0.10) this yields: n+1 


s + y/s 2 -J- 1 




s 2 + cos 2 - 

1 n+1 


(3.0.14) 


(_l)j+i . n j T j x 

— sm sin — — 

n + 1 n + 1 n + 1 


sin 2 - 2 ~- 

n+ 1 

n+1 


For (3.0.13) we note that cos 
ail but two terms in the sum: 


— cos 


(n+1— j)x 
n+1 * 


Therefore, the limit process eliminates 


(3.0.15) 


1 

n + 1 



jrx 

n + 1 


+ - sin 2 

Li 


(n + 1 — j)x \ 
n + 1 }' 


Expressions (3.0.14) and (3.0.15)are clearly equal. This, combined with the fact that the 
numerators have degree n - 1, implies that the two functions are the same.o 

Using the error formula for Gaussian quadrature we directly obtain an error estimate 
for the approximation, K n (t ), to S(t): 
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Theorem 3. For 0 < T < oo: 


( /y\2»+i 11 \ 

j-—, -r+||5|| £i(10i0o)) ). 


Proof: Noting that: 


(3.0.16) 
we have: 


M*)i < 


1 

n + l 


n 


T. 

i=i 


lx 

n + 1 



— w 2 dw = 


1 

2’ 


(3.0.17) \(Mt)/t) - K n {t)\ < min (j Q)’” l + |(7 l( t) /t) |j . 

Integrating this expression yields the desired result.o We note that this result clearly implies 
the rapid convergence of the approximation for fixed T as n is increased. 

This result can be extended to approximation of the pseudodifferential operator a 1 / 2 
in higher space dimensions. Introduce the spaces L 2 (H p , [0,T]) with norm: 

T 

(3 '°- 18) IKm = l f R Ji + |*| 2 F|«(*, t)\*dkdt. 

Writing: 

(3 - 0 ' 19) al/I = §i+S(\k\,t), K-=^+K n (\k\,t), 

we have: 


Theorem 4. For any e, l > 0, 0 < T < oo, there exists N such that for any n > IV 
and u 6 L 2 {H p , [0, T]); y 

||(a 1/2 - B n )u\\ p _ 2 _ g [0 T] < e|[u|J p ^ 0 >r j. 


Proof We first estimate the error in Fourier space. Note that: 

(m) *•<*■*> ■ S g*- ■ - ((«-ipi.). 


Then: 

(3.0.21) 


ll(<£ K-n) * u||Lj([0,r)) 5 II* - ^n||£ I ([0,TJ)l|«||£ 3 ([0,r]). 


The first term on the right can be estimated using the result of Theorem 3 and the 

change of variables z = |Jb|t: 


il* 


Knlkao.rj) - 1*1 f* T \(Mz)/z) - g sin 2 cos ( (cos ldi 


(3.0.22) 


< min M*l 




m 2n+1 __JL.„ i,, 2 


V 2 / ( 2 n + 'l)!’ 2 ^ T + l fc lll( J i(’)/('))IUi([o l00 )) I - 
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Given e, 8 > 0 we can clearly choose n sufficiently large that: 

(3 -°‘ 23) (! + |4| 2 )-< 1 +<*/ 2 »|| < S - K„|| il(( „ ir)) < e . 

Finally: 

ll<? - K n \[ < j R J l + |4| 2 ) p (l + |*p)-(2+^| jc 5 _ *„||i l( , wl) ||fi||J i(M)( ft 
(3.0.24) < e 2 ||«|| 2 |[0|r j, 

completing the proof.o 

We are confident that this convergence result can be translated into a convergence 
result for the artificial boundary problem for smooth solutions of the wave equation in 
simple domains, using for example the stability results of Ha-Duong and Joly [9] This 
along with some numerical experiments, will be the sub ject of future work. 

4. Long Time Approximations. The error estimates discussed above are strongly 
^pendent °n T and break down as T - oo. For some applications, for example the so- 
e hnutmg amplitude problem [6, 11], or applications to viscoelasticity [141, lone time 
approximations are desired. Also, such approximations will have better large |Jb| behavior for 
midtidimensional problems. To achieve this we must approximate the convolution kernel in 
Ai([0, oo)) by functions in our class H. This greatly constrains the functions at our disposal 
as seen in the following elementary result: 

he U ” “ elem€nt ° /Xl([0 ’ OO)) ^ lfMp0hS ° }f 

Proof: Simply note that the only functions in K which axe in £,([0, oo)) take the form of 
(sums of) polynomials multiplying exponentials where the exponent has negative real part. 
Ihe transform of such a function has a pole in the left half plane.o 
An immediate corollary of this result is: 

Corollary 2. A function r e ft-fl^oo)) cannot have a transform which may be 
written as the ratio of polynomials of definite parity. 

For the multidimensional operator associated with i/s 2 + |4| 2 , our corollary implies that 
ong time approximate conditions, in the sense described here, cannot be local in space. For 
the balf-space problem they generally involve the operator whose symbol is 141. This is the 
so-called Dirichlet-to- Neumann map for the Laplace equation in the half space. The use 
of this operator for long time solutions of the wave equation has been suggested in [6] and 

[llj. A key practical issue is the efficient application of the nonlocal map. For progress in 
this regard see [8, 16]. * 6 

4.1. Laguerre Expansions. The kernels we wish to approximate, in particular «S(t), 
re generally m £ 2 ([0 ,qq)). In £ 2 , convergent approximating sequences are easily con- 
structed via expansions in orthogonal functions. A complete set of orthogonal functions 
whose elements are in 7Z are the Laguerre functions: 


C„(yt) = e %L n (7t), £ n 


_ (*-!)" 

(s + 2 )"+! 


(4.1.25) 
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m 

9m 

11 X] QmJ-'m 

1! YhQmJJ'm «5||i 2 (fo,oo)) 

0 

0.75 

0.87 

0.22 

1 

-0.15 

0.79 

0.19 

2 

-0.186 

0.60 

0.11 

3 

1.992 x 10- 2 

0.59 

0.11 

4 

9.4512 x 10“ 2 

0.48 

7.2 x 10- 2 

5 

2.08554 x 10~ 2 

0.48 

7.0 x 10“ 2 

6 

-4.94016 x 10~ 2 

0.42 

5.7 x 10-2 

7 

-3.34492 x 10- 2 

0.41 

5.0 x 10-2 

j 8 

1.96034 x 10~ 2 

0.38 

4.7 X 10-2 

Le 

3.24971 x 10“ 2 

0.35 

3.9 x 10-2 

10 

3.04379 x 10~ 4 

0.35 

3.9 x IQ- 2 


Table 1 


Coefficients and L p errors for Laguerre approximations to S , y — 3/2. 


Here, L n (z) is the nth Laguerre polynomial and 7 > 0 is a parameter which may be chosen 
to optimize the convergence rate. The expansion of a kernel Q is given by: 


(4.1.26) 


0 « X) 9n£n(lt), 

n — 0 


9n~l[ C n (^t)Q(t)dt. 

Jo 


We have written a program in Maple to compute these coefficients for arbitrary 7 and the 
kernel J 1 (t)/t. We have numerically computed the L\ errors in the resulting approximations. 
We find 7 = 3/2 to be a convenient choice, both from the point of view of convergence of 
the series and simplicity of expansion coefficients. The results are given in Table 1. For 
reference we note that the Ly norm of S is about 1.6 and the norm is about .65. 

It is evident that the convergence rate of this series in L\ is slow at best. Indeed, 
Laguerre series are not generally convergent in L\. (See [17].) We do not even have a proof 
in this particular case. We have not yet implemented boundary conditions based on these 
expansions, but plan to do so in the future. 


4.2. Exponential Expansions. Another convenient set of orthogonal functions in 1Z 
are exponential polynomials: 


(4-2.27) VJnt) = e-?P n ( 2e“* - 1), 

where P n (z) is the nth Legendre polynomial and 7 > 0 is again a free parameter. The 
orthogonality of these functions is easily established by mapping [0, 00) to [—1,1] by z — 
2e“" Y< — 1. We then expand a kernel Q by: 

m 

(4.2.28) Q « ]T g n Vn{lft), g n ~ 7(2 n + 1) f VJat)Q(t)dt. 

n=0 

The Laplace transform of V n takes the form: 


^ = E 




f-r s -1. Bffih ’ 

i=0 * r o 


(4.2.29) 
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m 

9m 

II T,9mVm ~ ffO.ooD 

II YlQmPm - «5||i 2 rro.oo'n 

0 

1.80998 x 10" 1 

2.1 

0.51 

1 

3.49884 x 10" 1 

1.2 

0.24 

2 

1.49008 x lO” 1 

0.81 

0.19 

3 

-1.46339 x 10- 1 

0.76 

0.14 

4 

-1.24158 x 10- 1 

0.59 

0.10 

1 & 

I 

9.35578 x 10" 2 

0.52 

8.1 x 10" 2 

# 

4.84350 x 10" 2 

0.50 

7.5 x lO’" 2 

7 

-8.60302 x 10~ 2 

0.43 

5.7 x 10“ 2 I 

8 

2.89748 x 10~ 2 

0.41 

5.4 x 10~ 2 

9 

3.50595 x 10~ 2 

0.41 

5.1 x 10~ 2 

10 

-5.67689 x lO” 2 

0.38 

4.3 x 10“ 2 


Table 2 “ 

Coefficients and L p errors for exponential Legendre approximations to S, 7 = 1 / 5 . 


where the coefficients a n j are easily tabulated. We have written a Maple .procedure to 
jcompute the expansion coefficients for S, and have tabulated the results in Table 2. We found 
that y — 1/5 was a reasonable choice from the point of view of speed of convergence. 

Again the convergence of the series is rather slow, particularly in L u and the errors are 
generally larger than for the corresponding term in the Laguerre series. Moreover, we have 
no convergence proof. 

4.3. Other Approximations. Given the disappointing convergence behavior of the 
series above, we are led to consider other means for constructing long time appro xima tions 
One possibility is to consider the integral representation of S used in our analysis of the 
Pade approximants (3.0.11). Setting 2 = iu and rewriting it as a contour integral along 
the imaginary axis, we can then deform the contour into the left half plane. A quadrature 
scheme will then produce approximations which decay exponentially in t. We plan to 
investigate this procedure in future work. 
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